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Abstract 

We describe the advantages and disadvantages of numerical methods when Bohmian trajectory- 
grids are used for numerical simulations of quantum dynamics. We focus on the crucial non- 
crossing property of Bohmian trajectories, which numerically must be paid careful attention 
to. Failure to do so causes instabilities or leads to false simulations. 



1 Introduction 

Bohmian Mechanics. Bohmian Mechanics [1] is a mechanical theory for the motion of particles. 
We describe the theory in units of mass m — h — 1 . Given the Schrodinger wave function ip t of an 
TV- particle system, the trajectories of the N particles g fc (g°, -0°; G R 3 , 1 < k < N, are solutions 
of 



dt 



= 3m 



: (9i,-,9jv) 



(1) 



with q k (t = 0) = ql as initial conditions and ip* denoting the complex conjugate of ip, so that 
ip*ip = V| 2 We denote the configuration point x = {x\, xjy) 6 R 3N , V x = (V Xl , ■ ■ • , V Xiv ) 
where \7 Xk is the gradient with respect to Xk £ M 3 . ip t is the solution of the Schrodinger equation 



. diptjx) 
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with the initial condition ipt=o(x) — ipt=o( x i> f n) = '0 O ( a; ii xjq). Using configuration space 
language the Bohmian trajectory of an N particle system is an integral curve q(t) £ R 3Ar of the 
following velocity field on configuration space 



v*(q,t)=Jm 



(V?V^t)(?) 



(^)(«) 



(3) 



dqjt) 
dt 



= f *(</,*) 
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Under general conditions one has global existence and uniqueness of Bohmian trajectories, i.e. the 
integral curves do not run into nodes of the wave functions and they cannot cross [5] . From now 
on we shall only talk about the trajectories as integral curves in configuration space. Note that for 
one particle the configuration space is equal to physical space, for more than one particle this is not 
the case. 
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1 Interpreting ip*ip as inner product, the generalization to spin wave functions is straight forward. 
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The empirical import of Bohmian Mechanics arises from equivariance of the \ip\ -measure: One 
readily sees that by virtue of ([3| the continuity equation for the Bohmian flow on configuration 
space is identically fulfilled by the density pt = \ipt\ 2 , then known as the quantum flux equation, 

-Q t hV x [\i/3{x,t)\ v w {x,t)) = 

This means that if the configuration of particles q = (q®, q^) is distributed according to 
Pa = ■■■Xn)\ 2 at time t = 0, then the configuration q(t) = {Qx(Qi,ip°',t), q N (q%,ip°;t)) 

is distributed according to p t = \ijj t (xi, ...ajjv)| 2 at any time t. Because of this, Bohmian Mechanics 
agrees with all predictions made by orthodox quantum mechanics whenever the latter are unam- 
biguous [H[3]. 



Hydrodynamic Formulation of Bohmian Mechanics. Write ip in the Euler form 

where R and S are given by real- valued functions. Equation ((3]) together with equation ((2]) separated 
in their real and complex parts gives the following set of differential equations 

1-v.sm 

^M = -Iv,(JJ(»,i)V,S( a; ,t)) 

This set of equations can be examined along the Bohmian trajectory q. We then obtain 

dq 



dt 



V q S(q,t) (5) 



^1 = -l R{q ,t)^S(q,t) (6) 
8S(q,t) l(dq\ 2 F(g) ( l V 2 q R(q,t) (?) 



dt 2 \dt J ^ J 2 R(q,t) 

Numerical Integration. Numerically this set of equations (JS])-© can readily be integrated [6] 
and offers some advantages over standard techniques for solving the Schrodinger equation (J2J) nu- 
merically. The basic idea, due to Wyatt [6], is that the Bohmian configuration space trajectories 
define a co-moving grid (Bohmian grid) in configuration space which is best adapted for the compu- 
tation of i/j. If initially the grid points (n points in M. 3N ) are \ip°\ 2 distributed, they will remain \ipt\ 2 
distributed for all times t. So the co-moving grid spreads dynamically according to the spreading of 
\"4>t\ 2 and the grid points will primarily remain in regions of space where |t/>t| 2 is large while avoiding 
regions of nodes or tails of |-0t| 2 which are numerically problematic. Therefore Wyatt's idea is this: 
Integrate the equations |S |l -([7 jl simultaneously, i.e. get the best adapted moving grid along with ifi, 
instead of integrating the Schrodinger equation on some fixed or otherwise determined grid. This 
is why the algorithm is particularly interesting for long-time simulations which usually demand 
huge computational effort on fixed grids. Therefore if one aims at the relevant parts of the wave 
function (where probabilities are high) one can use a fixed number n of \jj)\ 2 distributed Bohmian 
grid points in R 3N , so that the Bohmian grid simulation scales with the number N of particles [6] 
while conventional grid methods mostly scale exponentially with N. We say more in section [3] on 
how to distribute n grid points in a \ip\ 2 manner. 
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In order to perform the numerical integration of the set of differential equations ©-(J?]) we 
follow the straight-forward method described in [6]. The only crucial part in this is computing 
the derivatives involved in the set of differential equations. There are several techniques known 
and [6] gives a comprehensive overview. Among them one technique called least square fitting is 
commonly used. In this letter we argue that this method is inappropriate for integrating ©-([7]) for 
general initial conditions and potentials. Bad situations arise whenever Bohmian trajectories move 
towards each other, because least square fitting will allow crossings of the simulated trajectories 
which are not allowed for Bohmian trajectories. Once a crossing is encountered in a numerical 
simulation further computation can be aborted because this numerical wave function would differ 
immensely from the solution of the Schrodinger equation ((2). This will happen generically, i.e. 
for non-gaussian wave functions. To illustrate our argument we shall present a typical numerical 
example in one dimension with one particle. Therefore in what follows N = 1 and the n trajectories 
which we shall consider (making up the grid) are the possible trajectories of this one particle only. 



2 Least Square versus Polynomial Fitting 

In order to compute the derivatives encountered in ©-(J?]) we only consider two different types 
of fitting, the least square fitting and the polynomial fitting. Most other fitting algorithms are 
descendants of one or the other. Both provide an algorithm for finding e.g. a polynomial of 
degree, say (m — I), which is in some sense to be specified close to a given function / : K — > R 
known only on a set of, say n, pairwise distinct data points, (xj, /(xj))i<i< n , as subset of the graph 
of /. The derivative can then be computed from the fitting polynomial by algebraic means. Let for 
the further discussion 

V ■= {f(xi))x<i<n 

X : = (.Tp 1 ^)l<i<n,l<i<m 
a = ( a i)l<i<m £ K" 1 

S = (5i)i<i<„ G K" 

Using this notation the problem of finding a fitting polynomial to / on the basis of n pairwise 
distinct data points of the graph of / reduces to finding the coefficients of the vector a obeying the 
equation 

y = X ■ a + 5 

such that the error term 5 is in some sense small. Here the dot • denotes matrix multiplication. 
The fitting polynomial is given by p(x) = X)J=i OLjX 1 . 

Polynomial Fitting. For the case n = m we can choose the error term 5 to be identical zero 
since X is an invertible square matrix and a — X^ 1 ■ y can be straight-forwardly computed. 



Least Square Fitting. For arbitrary n > m there is no unique solution anymore and one needs 
a new criterion for finding a unique vector a. The algorithm of least square fitting uses therefore 
the minimum value of the accumulated error A := w(x — Xi)Sf where w : R — > K specifies a 

weight dependent on the distance between the i-th data point Xi and some point x where the fitting 
polynomial shall be evaluated. The vector a can now be determined by minimizing A as a function 
of a by solving 

= ( -2j2w{x - Xi)(y - X ■ a^r 1 j =0 

1<?'<™ V i=l / l<j<m 

2 Or more general an element of a m dimensional vector space of functions. The coefficients of the design matrix 
X determine which basis functions are used. 
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Note that for m = n and any non-zero weight w the algorithm of least square fitting produces the 
same a as the polynomial fitting would do since for the a determined by polynomial fitting the 
non-negative function A (a) is zero and thus a naturally minimizes the error term. 

Why Least Square Fitting is inappropriate for Bohmian Grids. The algorithm of least 
square fitting is well known for its tendency to stabilize numerical simulations by averaging out 
numerical errors. But exactly this averaging makes it hard to keep the grid points from crossing 
each other. In order to understand what happens during a numerical simulation recall the form of 
the equations of motion I©-® which have to be integrated step by step. The change in time of 
the phase S is determined by three terms in equation of 0. The first two terms form the classical 
Lagrangian and depend on the current velocities of the grid points and on the potential V. The 
third term is the so-called quantum potential and depends only on R and its second derivative. The 
quantum potential is what prevents Bohmian trajectories from crossing each other. It therefore 
needs special attention during the numerical integration. Now imagine initial conditions such that 
two grid points approach each other. An increase of the density of the grid points in a region where 
the two grid points move towards each other will cause a small bump in R. By small bump we mean 
a bump of the wave function shape on a microscopic scale, i.e. a scale defined by few grid points. 
Such a small bump of R may have large derivatives changing the quantum potential in ||7J). Thus, 
it is absolutely vital for a numerical simulation to implement a fitting algorithm that reconstructs 
not only R, respectively S, in an accurate way but also its second derivative. The tendency of the 
least square fitting algorithm is to average out those small bumps in R, respectively S. Hence, the 
simulation is blind to recognize grid points moving towards each other and hence does not prevent 
a crossing of these trajectories. 

Note that since the averaging occurs on a microscopic scale, situations in which grid points move 
only very slowly or do generically not move towards each otheJl are often numerically doable with 
least square fitting. On the other hand numerical simulations with least square fitting in general 
situations like the one discussed above are bound to break down as soon as two grid points get too 
close to each other. 

Why Polynomial Fitting is more appropriate for Bohmian Grids. We stressed above 
that small bumps may exponentiate numerical instability. A good numerical method must take 
note of the small bumps and must prevent their increase. Such a numerical method is provided by 
polynomial fitting, since there the polynomials go through all grid points. Therefore polynomial 
fitting recognizes the bumps of R and/or S and hence the resulting quantum potential recognizes 
the approaching grid points. Now recall the physics of the Bohmian evolution, which as we stressed 
in the introduction prevents trajectories from crossing each other. Therefore we expect that this 
method is self-correcting and hence stabilizing. 

The Boundary Problem. It has to be remarked that polynomial fitting creates a more severe 
problem at the boundary of the supporting grid than least square fitting. This problem is of 
conceptual kind since at the boundary there is a generic lack of knowledge of how the derivatives 
of R or S behave, see figured! 

We describe briefly how this conceptual problem can cause severe numerical instability. For this 
suppose that only the last grid point in the upper left plot of figure [1] is lifted a little bit upwards by 
e.g. some numerical error. Then the resulting quantum potential will cause the last grid point to 
move towards the second last one. This increases the density of grid points and thus again increases 
R in the next step such that this effect is self- amplifying (in fact growing super exponentially) and 
will effect the whole wave function quickly. 

At the boundary the good property of polynomial fitting, namely to recognize all small bumps 
works against Bohmian grids techniques. In contrary least square fitting simply averages these small 

3 For example a free Gaussian wave packet or one approaching a potential creating only soft reflections. 
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numerical errors out, see lower right plot in figure [TJ The conceptual boundary problem, however, 
remains and will show up eventually also with least square fitting. 

The Numerical Simulation. To demonstrate our argument we shall now give a numerical ex- 
ample for one particle in one dimension. We remark as mentioned earlier that Gaussian wave 
packets are unfit for probing the quality of the numerical simulations. Therefore we take as initial 
conditions two superposed free Gaussian wave packets with a small displacement together with a 
velocity field identically to zero and focus the attention on the region where their tails meet, i.e. 
where the Bohmian trajectories will move towards each other according to the spreading of the 
wave packets, see figure^ The physical units given in the figures and the following discussion refer 
to a wave packet on a length scale of lA = 10 -10 m and a Bohmian particle with the mass of an 
electron. 

The numerical simulation is implemented as in [6] with minor changes. It is written in MATLAB 
using IEEE Standard 754 double precision. In order to ease the boundary problem we use a very 
strong least square fitting at the boundary. Note that this only eases the boundary problem and 
is by no means a proper cure. In between the boundary, however, the number of data points for 
the fitting algorithm and the degree of the fitting polynomial can be chosen in each run. In this 
way it is possible to have either polynomial fitting or least square fitting in between the boundary. 
The simulation is run twice. First with the polynomial fitting and then with the least square fitting 
algorithm in between the boundaries. 

We take 7 basis functions for the fitting polynomial in both cases. Note that the choice should 
at least be greater or equal to 4 in order to have enough information about the third derivative of 
the fitting polynomial. In the first run with polynomial fitting the number of grid points used for 
fitting is 7 and in the second run for the least square fitting we choose 9 which induces a mild least 
square behavior. For the numerical integration we have used a time step of 10 _2 fs with the total 
number of 51 grid points supporting the initial wave function. Please find the source code of our 
numerical simulation at the end of this letter. 

Results. The first run with polynomial fitting yields accurate results and does not allow for 
trajectory crossing way beyond 5000 integration steps, i.e. 50fs, see figures [3] and |4j The second 
run with least square fitting reports a crossing of trajectories already after 430 integration steps, 
i.e. 4.3fs, and aborts, see figures [5] and [6] The numerical instability can already be observed earlier, 
see left-hand side of figure An adjustment of the time step of the numerical integration does not 
lead to better results in the second run. The crossing of the trajectories of course occurs exactly in 
the region were the grid points move fastest towards each other. Referring to the discussion before, 
the small bumps in R created by the increase of the grid point density in this region is not seen 
by the least square algorithm and thus not seen by the numerical integration of equations ©-([7]), 
compare figure [3] and [5j The approaching grid points are not decelerated by the quantum potential 
and finally cross each other, see figure while in the first run they begin to decelerate and turn 
to the opposite direction during the time between 3fs to 5fs, see figure [3J To visualize how the two 
algorithms "see" R 2 and the velocity field during numerical integration these entities have been 
plotted in figures [3] and [5] by merging all fitting polynomials in the neighborhood of every grid point 
together (from halfway to the left neighboring grid point to halfway to the right neighboring grid 
point). In figure [5] one clearly observes the failure of the least square fitting algorithm to see the 
small bumps in R. 

Note that in some special situations in which the time step of the numerical integration is chosen 
to be large, polynomial fitting may lead to trajectory crossing as well. This may happen when the 
number of time steps in which the simulation has to decelerate two fast approaching grid points is 
not sufficient. This effect is entirely due to the fact that numerical integration coarse grains the 
time. The choice of a smaller time step for the simulation will always cure the problem as long as 
other numerical errors do not accumulate too much. 

The relevant measure of quality of the numerical simulation is naturally the L 2 distance between 
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the simulated wavefunction Re tS and the analytic solution of the Schrodinger equation %jj t , i.e. 

1 /2 

(/ dx \ipt{x) — R(x, t)e lS ^ x '^\ 2 ) , and is spelled out on the left-hand side of the figures Q] and [5] 
for both runs. 

3 Using \i/j°\ 2 as Initial Distribution. 

As discussed in the introduction, if the grid point are distributed according to IV' ! 2 , they will 
remain so for all times. Doing so will bring two advantages but comes at the price of a more 
severe conceptual boundary problem as discussed in the last section. The first advantage is that 
this Bohmian grid avoids regions where R is very small and thus where the computation of the 
quantum potential becomes numerically tricky. The second advantage is that now the system is 
over-determined and the actual density of grid points must coincide with R 2 for all times. This 
could be used as an on the fly check whether the grid points behave like Bohmian trajectories or not. 
If not the numerical integration does not give a good approximation to the solution of ((5])- (0). It 
could also be considered to stabilize the numerical simulation with a feedback mechanism balancing 
R 2 and the distribution of the particle positions which may correct numerical errors of one or the 
other on the fly. 

One way of choosing \ip°\ 2 distributed initial positions for the grid points q® for j = 1 . . .n for 
some large number n is the following. Choose q® in such a way that 

^o (g o )|a g? + i-g?-i = I (8) 

This formula can be used to compute the grid points iterativeljQ starting with some grid point near 
the maximum of \tp(q, 0)| 2 . Another way is to simulate I?/" ! 2 distributed random variables q® by 
the commonly used techniques. Here, however, it has to be taken care that the randomly chosen 
g° do not lie too close together. If they do, either some of them must be deleted or the time step 
of the numerical integration must be adjusted carefully. A rule of thumb is that the distance of 
approaching grid points divided by their relative velocity has to be a lot smaller than the chosen 
time step. 

4 Conclusion 

Using polynomial fitting instead of least square fitting in between the boundary increases the 
stability of the numerical integration immensely. This is due to the fact that polynomial fitting in 
contrary to least square fitting does not average over the microscopic structure of the function to 
fit and therefore reconstructs needed derivatives more accurately. The discussed boundary problem 
arises more severely for polynomial fitting because of this fact. However, this problem is generic 
to the numerical integration considered here and also occurs using least square fitting but in a 
milder way. We suggest that the smoothness of the decay of the wave function near the boundary 
should be the guide for solving this problem which requires a detailed study. Furthermore, we have 
discussed that the Bohmian grid is best adapted to the problem of numerical integration of ||5j)-|(7J) 
because the grid points naturally avoid regions where R becomes very small. 
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Figure 1: Fitting polynomial of sixth degree and its second derivative near the boundary. Upper left: 
the function to fit smoothly decays. Both least square and polynomial fitting give identical fitting 
polynomials mimicking this decay. Upper right: polynomial fitting was used while the function 
value at fourth grid point was shifted upwards by I 0~ 4 (not visible in plot) to simulate a numerical 
error in order to see how sensitive polynomial fitting reacts to such an error. Lower left: all values 
at the grid points were shifted by an individual random amount of the order of 1CP 3 and polynomial 
fitting was used. Lower right: all values at the grid points were shifted by the same random numbers 
as on the lower left while the additional grid points were also shifted by individual random amounts 
of the order of f0~ 3 and least square fitting was used. By comparison with upper left plot one 
observes the robustness of least square fitting to such numerical errors at the boundary. 
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Figure 2: Left: Initial wave function at time zero. Right: Its velocity field after a very short time. 
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Figure 3: Simulated |V>| 2 together with the velocity field at times t = 3.8fs and t = 15fs using 
polynomial fitting in between the boundary. Left: the center part of the wave function where grid 
points move towards each other. Right: the whole wave function at a later time. The grid points 
have all turned and move apart from each other. The kinks at ±20 • 10~ 15 to in the velocity field in 
the lower right plot are due to the transition from least square fitting at the boundary to polynomial 
fitting in between the boundary (see paragraph: The Boundary Problem in section [2|). 




Figure 4: Both plots belong to the simulation using polynomial fitting in between the boundary. 
Left: L 2 distance between the simulated wavefunction and the analytic solution of the Schrodinger 
equation. Right: A plot of the trajectories of the grid points. Note how some trajectories initially 
move towards each other, decelerate, and finally move apart. 
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Figure 5: Left: same as left-hand side of [3] but this time using least square fitting throughout. 
Right: a short time later. Note in the upper left plot, the fitting polynomials of least square fitting 
fail to recognize the relatively big ordinate change of the grid points. So the grid points moving 
towards each other are not decelerated by the quantum potential. Finally after t = 4.3fs a crossing 
of the trajectory occurs. 
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Figure 6: Both plots belong to the simulation using least square fitting throughout. Left: L 2 
distance between the simulated wavefunction and the analytic solution of the Schrodinger equation. 
Right: A plot of the trajectories of the grid points. Note by comparison with right-hand side of 
figure [4] the failure of least square fitting to recognize grid point trajectories moving towards each 
other until they finally cross and the numerical simulation aborts. 
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SOURCE CODE OF THE NUMERICAL SIMULATION: 



7,7, definition of constants and rendering of initial value 

7oCONSTANTS & SYMPBOLS***************************************************** 



f ovx 

numParticles 
dt 

numSteps 

stencil 

degree 



8; 
51; 

10~-2; 
6000; 

3; 

6; 



s_stencil 



s_degree 



= 3; 



= 6; 



°/»field of view on the real line [-fovx, fovx] 

°/ number of Bohmian particles. 

°/ n time step for numerical integration. 

°/ number of iterations of the numerical 

7,integration. 

°/ number of points left and right w.r.t. some 
Zpoint for fitting 0.5*log(rho) . 
"/.degree of polynomial for fitting 
Zplease note: stencil=degree/2 means a 
Zprecise polynomial fit of 'degree' degree 
%while stencil>degree/2 means a least square 
7,fit with a polynomial of degree 'degree'. 
°/»log(sqrt (rho) ) number of points left and 
Zright w.r.t. some point for fitting the 
Zphase . 

°/ degree of polynomial for fitting the phase. 
%the above note hols for the phase. 



sym_x = sym( ' sym_x' , ' real ' ) ; "/.dummy symbol for symbolic math 

% 

7 Analytic solution of an initial gaussian evolved according to the free 
7 Schroedinger equation. 

7, INPUT : t=time of evolution, x=position, sigma=initial width of the 
7o gaussian 

7o OUTPUT: complex value of the wavef unction 

gauss = @(t, x, sigma) (sigma/ (pi* (sigma+i*t) ~2) ) . . . 

. "(1/4) . *exp(-x. ~2/(2*(sigma~2+t~2))*(sigma-i*t)) ; 

% 

7, psi=analytic solution of a superposition of initial gaussians evolved 
7, according to the free Schroedinger equation. 
7, rho,S,vel=|psi| ~2, the phase, the velocity field of psi 
7« INPUT : t=time of evolution, x=position 
7« OUTPUT: real resp. complex function value 
psi = @(t, x) l/sqrt(2)*(gauss(t,x-3,4)+gauss(t,x+3,4)) ; 
rho = @(t, x) psi (t ,x) . *conj (psi (t ,x) ) ; 
S = @(t, x) imag(log(psi (t ,x) . /abs (psi (t ,x) ) ) ) ; 

vel = @(t, x) imag(subs(dif f (psi (t , sym_x) ,sym_x) . . . 
. /psi (t , sym_x) ,sym_x,x)) ; 

7.INITIAL DATA 



clear x; 



7«initial distribution of the 
7 Bohmian particles 
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x = [-f ovx : (2*f ovx/mimParticles) : f ovx] ; "/generate a uniform distribution 

"/plot the initial data: blue=initial rho and green=initial veleocity field 
plot(x, vel(0,x), 'g', x, rho(0,x), 'b.-'); 
axis ([-f ovx fovx -1 1]); 
grid on; 

77 bohmian propagation 
warning off ; 

"/Unitialize the array that keeps track of the (x,0.5*log(rho) ,S,time) 
"/data history. 

history = [{x vel(0,x) 0.5*log(rho(0,x)) S(0,x) 0}]; 

"/.loop over the numer of time steps of numerical integration 
for (step = l:numSteps) 

xlen = length (history{step, 1}) ; 

x_list = history{step, 1}; 

vel_list = history{step,2}; 

logsqrtRho_list = history{step,3}; 
s_list = history{step,4}; 

"/.first loop over the all Bohmian particles updating the 
°/x_list and s_list 
for (c = l:xlen) 

"/find stencil points for logsqrtRho 

left_point = c - stencil; 

right_point = c + stencil; 

deg = degree; 

if (left_point < 1) 
left_point = 1; 

right_point = stencil*2+l + round (numParticles/7) ; 
"/decrease degree at the boundary as dodgy fix to avoid 
"/bouadary problems 
deg = 2; 
elseif (right_point > xlen) 

left_point = xlen - (2*stencil+l) - round (numParticles/7) ; 
right_point = xlen; 

"/decrease degree at the boundary as dodgy fix to avoid 
"/bouadary problems 
deg = 2; 

end 

"/fit logsqrtRho_list 

fitdata_x = history{step, 1} (lef t_point : right_point) ; 
fitdata_y = history{step,3}(left_point :right_point) ; 
logsqrtRho_f it = polyf it (f itdata_x, fitdata_y, deg); 
"/compute derivatives 

d_logsqrtRho_f it = polyder (logsqrtRho_f it) ; 
dd_logsqrtRho_f it = polyder (d_logsqrtRho_f it) ; 
"/compute quantum potential 
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Q = -l/2*(polyval(dd_logsqrtRho_f it ,x_list (c) ) + ... 

polyval(d_logsqrtRho_f it ,x_list (c) ) . ~2) ; 
"/.update phase: S(t+dt)=S(t)+dt(l/2v(t)~2-Q(t)) 
s_list(c) = s_list(c) + (l/2*vel_list(c) .'2 - Q)*dt; 
%update position: r(t+dt)=r(t)+v(t)*dt 
x_list(c) = x_list(c) + vel_list(c)*dt; 

end 

7 second loop over the all Bohmian particles updating 
"/.the vel_list and logsqrtRho_list 
for (c = l:xlen) 

"/find stecil points for S 

s_left_point = c-s_stencil; 

s_right_point = c+s_stencil; 

deg = s_degree; 

if (s_lef t_point < 1) 
s_left_point = 1; 

s_right_point = s_stencil*2+l + round (numParticles/7) ; 
"/.decrease degree at the boundary as dodgy fix to avoid 
"/.bouadary problems 
deg = 2; 
elseif (s_right_point > xlen) 

s_left_point = xlen - (2*s_stencil+l) - round (numParticles/7) ; 
s_right_point = xlen; 

"/decrease degree at the boundary as dodgy fix to avoid 
"/bouadary problems 
deg = 2; 

end 

"/.fit s_list 

fitdata_x = x_list (s_lef t_point : s_right_point) ; 
fitdata_y = s_list (s_lef t_point : s_right_point) ; 
s_fit = polyf it (f itdata_x, fitdata_y, deg); 
"/.compute derivatives 
d_s_fit = polyder (s_f it) ; 
dd_s_f it = polyder (d_s_f it) ; 
°/.v(t+dt) = S' (t+dt) 

vel_list(c) = polyval(d_s_f it ,x_list (c) ) ; 
°/.R(t+dt) = R(t)*exp(-l/2 S"(t+dt)*dt) 
logsqrtRho_list (c) = logsqrtRho_list (c) - . . . 

l/2*polyval (dd_s_f it , x_list (c) ) *dt ; 

end 

"/compile the tuple (x,vel,0.5*log(rho) ,time) 
newSnapshot = [{x_list vel_list . . . 

logsqrtRho_list s_list history{step,5}+dt>] ; 
"/and save the snapshot 
history = [history; newSnapshot] ; 

"/.compute minimal distance of the Bohmian particles 

dxmin = min(x_list (2 : length(x_list) ) -x_list (1 : length(x_list) -1) ) ; 

if (dxmin <= 0) 

sprintf ('TRAJECTORY CROSSING OCCURED AT TIME: °/.f ' , . . . 
history{step , 5}+dt) 
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break; 

end 

if (step == 1 I mod(step/numSteps*100,10) == 0) 
%plot integrated velocity field 

plot (x_list (1) : . 1 : x_list (xlen) , vel(history{step,5}+dt , 
x_list(l) :0.1:x_list(xlen)) , >r' ,x_list ,vel_list , 'b 
axis ([-25 25 -5 5] ) ; 
grid on; 

sprintf ( 'Progress : '/,.lf%'/ , timestep: °/ f fs, min deltax: 

step/numSteps*100, history{step,5}+dt, dxmin) 
pause(0 . 1) ; 

end 

end 

°/o°/o render rho movie (red is the analytic solution) 

frame = 1; 

clear ( ' rhoMov ' ) ; 

for step=l : 10 : length (history) 

plot (-25:0. 1:25, rho(history{step, 5}, -25:0. 1:25) ,'r' ,. . . 
history{step , l},exp(history{step,3}*2) , 'b. ' , . . . 
history{step, 1},0, 'b*') ; 

axis ([-25 25 0.3]) ; 

grid on; 

pause (0 . 05) ; 

rhoMov(f rame) = getframe; 
frame = frame + 1; 

end 

°/o°/o render vel movie (red is the analytic solution) 

frame = 1; 

clear ( ' velMov' ) ; 

for step=l : 10 : length (history) 

plot (-20:0. 1:20, vel(history{step, 5}, -20:0. 1:20) , 'r' , . . . 
history{step, l},history{step,2>, 'b. ') ; 

axis ([-25 25 -5 5] ) ; 

grid on; 

velMov(f rame) = getframe; 
frame = frame + 1; 

end 

°/ °/o render velocity field seen by the above algorithm 
diagAxis = [-10 10 -1 1] ; 

°/«loop over the number of time steps of numerical integration 
for (step = l:length(history)) 

xlen = length (history{step, 1}) ; 

xvals = [] ; 
yvals = [] ; 
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'/.loop over all particles 
for (c = l:xlen) 

'/.find stecil points for S 

s_lef t_point = c-s_stencil; 

s_right_point = c+s_stencil; 

deg = s_degree; 

if (s_lef t_point < 1) 
s_left_point = 1; 

s_right_point = s_stencil*2+l + round (numParticles/7) ; 
"/.decrease degree at the boundary as dodgy fix to avoid 
'/.bouadary problems 
deg = 2; 
elseif (s_right_point > xlen) 

s_left_point = xlen - (2*s_stencil+l) - round (numParticles/7) ; 
s_right_point = xlen; 

'/.decrease degree at the boundary as dodgy fix to avoid 
'/.bouadary problems 
deg = 2; 

end 

'/.fit s_list 

fitdata_x = history{step, 1} (s_lef t_point : s_right_point) ; 
fitdata_y = history{step,4} (s_lef t_point : s_right_point) ; 
s_fit = polyf it (f itdata_x, fitdata_y, deg); 
'/.compute velocity 
d_s_fit = polyder (s_f it) ; 
'/.compute graph 

if (c == 1) 

xl = history{step,l}(c)-l; 
x2 = history{step, l}(c) ... 

+ (history{step, 1} (c+1) -history{step, l}(c) ) /2 ; 
elseif (c == xlen) 

xl = history{step, l}(c-l) ... 

+ (history{step, 1} (c) -history{step, 1} (c-1) ) /2 ; 
x2 = history{step, l}(c)+l ; 

else 

xl = history{step, l}(c-l) ... 

+ (history{step , 1} (c) -history{step , 1} (c-1) ) /2 ; 
x2 = history{step, l}(c) ... 

+ (history{step , 1} (c+1) -history{step, 1} (c) ) II ; 

end 

xvals = [ xvals xl : ( (x2-xl) /100) :x2 ]; 

yvals = [ yvals polyval (d_s_f it ,xl : ( (x2-xl) /100) :x2)] ; 

end 

plot (xvals, vel(history{step, 5}, xvals) ,'z' , xvals , yvals , 'b' , . . . 

history{step , 1} , history{step , 2} , ' b . ' ) ; 
axis (diagAxis) ; 
grid on; 
pause (0. 1) ; 

end 



15 



7,7o render rho field seen by the above algorithm 
diagAxis = [-30 30 0.145]; 

"/.loop over the number of time steps of numerical integration 
for (step = l:length(history)) 

xlen = length(history{step, 1}) ; 

xvals = [] ; 
yvals = [] ; 

'/.loop over all particles 
for (c = l:xlen) 

'/.find stencil points for logsqrtRho 

left_point = c - stencil; 

right_point = c + stencil; 

deg = degree; 

if (left_point < 1) 
left_point = 1; 

right_point = stencil*2+l + round (numParticles/7) ; 
"/.decrease degree at the boundary as dodgy fix to avoid 
'/.bouadary problems 
deg = 2; 
elseif (right_point > xlen) 

left_point = xlen - (2*stencil+l) - round (numParticles/7) ; 
right_point = xlen; 

'/.decrease degree at the boundary as dodgy fix to avoid 
'/.bouadary problems 
deg = 2; 

end 

'/.fit logsqrtRho_list 

fitdata_x = history{step, 1} (lef t_point : right_point) ; 
fitdata_y = history{step,3}(left_point :right_point) ; 
logsqrtRho_f it = polyf it (f itdata_x, fitdata_y, deg); 
'/.compute derivatives 

d_logsqrtRho_f it = polyder (logsqrtRho_f it) ; 
dd_logsqrtRho_f it = polyder (d_logsqrtRho_f it) ; 
'/.compute quantum potential 

Q = -l/2*(polyval(dd_logsqrtRho_f it ,x_list (c) ) + ... 

polyval(d_logsqrtRho_f it ,x_list (c) ) . ~2) ; 
'/.compute graph 

if (c == 1) 

xl = history{step, l}(c) -1 ; 
x2 = history{step, l}(c) ... 

+ (history{step , 1} (c+1) -history{step , 1} (c) ) /2 ; 
elseif (c == xlen) 

xl = history{step, l}(c-l) ... 

+ (history{step , 1} (c) -history{step , 1} (c-1) ) /2 ; 
x2 = history{step, l}(c)+l ; 

else 

xl = history{step, l}(c-l) ... 
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+ (history{step , 1} (c) -history{step , 1} (c-1) ) /2 ; 
x2 = history{step, l}(c) ... 

+ (history{step , 1} (c+1) -history{step, 1} (c) ) /2 ; 

end 

xvals = [ xvals xl:0.05:x2 ]; 

yvals = [ yvals polyval (logsqrtRho_f it ,xl : . 05 : x2)] ; 

end 

plot (xvals, rho(history{step, 5}, xvals) , 'r' , xvals, exp(2*yvals) , 'b' , . . . 

history{step, l},exp(2*history{step,3}) , 'b. ' ) ; 
axis (diagAxis) ; 
grid on; 
pause (0. 1) ; 

end 

7,7, render trajectories 

diagAxis = [-12 12 length (history) *dt] ; 
hold off; 
grid off; 

"/.loop over the grid points 
for (j = 1 :numParticles) 

xvals = [] ; 

yvals = [] ; 

for (step = 1 : length (history) ) 

xvals = [xvals; history{step, l}(j)] ; 
yvals = [yvals; history{step,5>] ; 

end 

plot (xvals , yvals , ' k- ' ) ; 
axis (diagAxis) ; 
hold on; 
pause (0. 1) ; 

end 

hold off; 

7,7, render L~2 error 
psi2error = [] ; 

for (step = l:length(history)) 
xdif f = [] ; 
ydif f 2 = [] ; 

for (j = 1 :numParticles-l) 

xdiff = [xdiff; history{step, 1} ( j+l)-history{step , 1} (j )] ; 
ydiff2 = [ydiff2; abs(exp(complex(0,history{step,4}(j+l))) . . . 

. *exp(history{step,3}(j+l) ) . . . 

- psi((step-l)*dt,history{step,l}(j+l)))~2] ; 

end 

psi2error = [psi2error; sqrt (sum(xdif f .*ydiff2))] ; 

end 

plot ( (1 : length (history) ) *dt ,psi2error , 'k' ) ; 
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